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ABSTRACT 

Cosmic microwave background (CMB) temperature anisotropies follow a Gaussian 
statistical distribution in the standard inflationary model, but there are non-Gaussian 
contributions due to astrophysical foregrounds. The detection of the non-Gaussian 
signal due to extragalactic point sources and its distinction from the possible intrinsic 
non-Gaussianity is an issue of great importance in CMB analysis. The Mexican Hat 
Wavelet Family (MHWF), which has been proved very useful for the detection of 
extragalactic point sources, is applied here to the study of non-Gaussianity due to point 
sources in CMB maps. We carry out simulations of CMB maps with the characteristics 
of the forthcoming Planck mission at 70 and 100 GHz and filter them with the MHWF. 
By comparing the skewness and kurtosis of simulated maps with and without point 
sources, we are able to detect clearly the non-Gaussian signal due to point sources for 
flux limits as low as 0.4 Jy (70 GHz) and 0.3 Jy (100 GHz). The second and third 
members of the MHWF perform better in this respect than the Mexican Hat Wavelet 
(MHW1) and much better than the Daubechies 4 wavelet. We have also estimated 
the third order, K^, and fourth order, K4, cumulants produced by point sources at 
these Planck channels by means of a fit with the MHWF. The average relative errors 
with respect to the real values are below 12% for fluxes down to 0.6 Jy (70 GHz) 
and 0.4 Jy (100 GHz). The values of these cumulants allow us to distinguish between 
different source counts models. From the estimated cumulants and assuming a power 
law for the source number counts we are able to obtain the coefficients a and A of 
the differential number counts, a = 2.19 ± 0.46 {a = 2.26 ± 0.19), A = 24.3 ± 4.3 
(A = 21.0 ± 4.2) for 70 (100) GHz, assuming a flux limit of 1 Jy. These results are 
consistent with the values obtained from simulations containing only point sources, 
which are: a = 2.32 ± 0.06 (a = 2.35 ± 0.16), A = 22.1 ± 1.5 {A = 20.0 ± 3.7) for 70 
(100) GHz. 
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1 INTRODUCTION 

According to the standard cosmological model, CMB tem- 
perature anisotropies are caused by the inhomogeneities 
in the distribution of matter and radiation present at 
the decoupling time. These inhomogeneities are the prod- 
uct of quantum fluctuations amplified in the inflation- 
ary era and follow a Gaussian statistical distribution. 
Therefore, the observed CMB anisotropies are a realiza- 
tion of a homogeneous and isotropic Gaussian random 
field on the sphere. Nevertheless, some degree of non- 
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Gaussianity ca n be introduced, for instance, by non stan- 
dard inflation jBartolo et al.| |2004l) or by t opological de- 
fects IITurok and Spe rgel 19901 iDurrerl 1999 h The analysis 
of the first and three-year Wilkinson Microwave Anisotropy 
Probe (WMAP) data by the WMAP team shows that 
CMB temperatur e fluctuations are consistent with a Gaus- 
sian distribution jKomatsu et alJl2003l : ISoergel et al.ll2006ft 
in agreement with the standard inflation paradigm. How- 
ever, some studies have detected non-Gaussian featur es in 
the WMAP data, for in stance in IVielva et all i2004ft and 
ICruz et alJ <2005l l200d) a cold non- Gaussian spot of un- 
known origin was detected and studied by using the Spher- 
ical Mexican Hat Wavelet. Non- Gaussian signatures have 
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also be en found by using different methods : phase corre- 
lations dChiang et alJl2003tJColes et al.ll2004). local curva- 
ture methods llHansen et all2004l:ICa bella et al.l2005l) , mul- 
tivariate analysis (Dinccn and Coles 2005), Minkowski func- 
tional <Park|2004|lEriksen et alJl2004) and higher criticism 
n et al.l l2(] 

On the other hand, astrophysical foregrounds con- 
taminate the cosmological signal and also give rise to a 
non-Gaussian contribution. In particular, the emission due 
to point sources (radio and infrared galaxies are seen as 
point-like objects due to their very small projected angu- 
lar size in comparison with the ty pical experimental res- 
olutions) is clearly non-Gaussian l|Ko-ma tsu and gper^el I 
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20oj iHerranz et al J 1200 i iBarreiro et al Jl2()()(f). "in the last 
few years, new high frequency (v ^ 10 GHz 



surveys 



of radio sources have been carried out dWald ram et alJ 
120031: iRicci et al.ll20oi) . All the recently published data on 
counts and on the emission properties of extragalactic ra- 
dio sources have guided the development of new cosmo- 



dicting number counts and source statistics at high radio fre- 
quencies llDe Zotti et alJl2005l) . We will use this new model 
to carry out point source simulations. 

It is very important to assess the degree of non- 
Gaussianity generated by point sources in order to distin- 
guish it from the non- Gaussian sig nal coming from other 
origins. The ESA Planck satellite (iMandolesi et alJ Il99st 
iPuget et aDll998l) . which will be launched in 2008, will an- 
alyze the CMB anisotropics with unprecedented accuracy 
and resolution and will make possible a very precise non- 
Gaussianity study, for which the use of wavelets can be very 
suitable. 

The Mexican Hat Wavelet (MH W1) has bee n very 
useful in detectin g point sources JCavon et al] 120001 
IVielva et all l2001allb|) in CMB maps, by using the signal 
amplification going from real space to wavelet space. The 
Mexican Hat Wavelet Family (MHWF) was introduced by 
iGonzalez-Nuevo et al] <|2006f) as an extension of the stan- 
dard MHW1, which produced higher amplification (see be- 
low) and then a better performance at point source detec- 
tion. These new wavelets we re also used for point sourc e 
detection in Planck maps in lLopez -Canieeo et al~l (2006). 
The MHWF is a family of isotropic 2D wavelets obtained 
by applying the Laplacian operator iteratively to the 2D 
Gaussian. 

The expression of these new wavelets in real space is 



(1) 
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where x is the radial coordinate and <p is the 2D Gaussian 



<p(x) = 



2tt 
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Note that ipi(x) is the standard MHW1, we call the other 
corresponding members of the family MHWn. Any member 
of the family can be written very easily in Fourier space as 
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IGonzalez-Nuevo et al] ]2006l) considered a point source con- 
volved with a Gaussian instrumental beam of dispersion 7 



7 = 



FWHM 
V 81o g( 2 ) 
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where FWHM denotes the Full Width Half Maximum of the 
instrumental beam. 

The point source was embedded in a CMB map, in- 
cluding noise and all the relevant foregrounds. Then, they 
filtered the map with the first members of the MHWF at 
different scales R. 



*„(x;6, R) = 



R 
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The amplification, A, is defined as the ratio between the 
intensity of the point source I w and the rms deviation a w 
in the wavelet filtered image divided by this ratio in the 
original image filtered with a Gaussian-shaped instrumental 
beam. 
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I w can be calculated for the different members of the MHWF 
as 
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with /3 = R/y. 

It was proved that when filtering CMB maps with the 
first three members of the MHWF at an optimal scale, R, the 
MHW2 and MHW3 produced a higher amplification than 
the standard MHW1 and were hence more effective at source 
detection. Other higher n members of the MHWF did not 
yield higher amplifications. 

After removing the detected point sources, there could 
still be a remaining non-Gaussianity due to undetected point 
sources. It is of great interest to characterize this non- 
Gaussianity in order to distinguish it from the intrinsic one. 

The detection of the non- Gaussian signal due to point 
sources poses a similar problem to that of point source de- 
tection, these new wavelets will enhance the signal (non- 
Gaussian) due to point sources, thus making it possible to 
detect this kind of non-Gaussianity more effectively. Since 
the skewness and kurtosis of a Gaussian distribution are zero 
and three respectively, any significant deviation from these 
values will mean the presence of a non-Gaussian signal, so 
that the skewness and the kurtosis of the wavelet filtered 
maps can be used as estimators of the non-Gaussianity due 
to point sources. 

Filtering CMB maps with the MHWF at optimal scales 
amplifies point sources and can also be used to determine 
the third and fourth order cumulants produced by the point 
sources. From these cumulants, we can also determine the 
parameters A and a, if we assume that the differential num- 
be r counts can be expressed as a power law dn/ds = AS~ a . 
In lPierpaolil i2003l) a chi-squared fit to a power law was ex- 
plored considering the moments of the CMB maps. We have 
considered in this paper a chi-squared fit of the third or- 
der and fourth order cumulants of the wavelet filtered maps 
combining several scales including the optimal ones. 

The plan of the paper is as follows: in §2, we calcu- 
late the cumulants of a wavelet filtered map from those of 
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the original map, this allows us to relate the skewness and 
kurtosis of the filtered maps to the amplification. In §3, we 
carry out realistic CMB simulations with and without point 
sources, considering the characteristics of the Planck exper- 
iment at 70 and 100 GHz and perform the non-Gaussianity 
detection by filtering with the MHWF. In §4, we obtain the 
third and fourth order cumulants produced by point sources 
by means of a fit of the cumulants of wavelet filtered CMB 
maps and we also determine the parameters a and A. Finally 
in §5, we draw the main conclusions of the paper. 



2 SKEWNESS AND KURTOSIS OF A 
WAVELET FILTERED MAP 

Skewness and kurtosis provide us with suitable tests for de- 
termining whether certain data come from a Gaussian dis- 
tribution or not. The skewness, s, is defined as 



Its 



(8) 



where /i3 is the third central moment of the distribution 
and a its rms deviation, the skewness is zero for a Gaussian 
distribution. The kurtosis, k, is the fourth central moment 
divided by the squared variance 



k- 
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and its value is 3 in the Gaussian case. We are interested 
now in calculating the skewness and kurtosis due to point 
sources in a wavelet filtered map from these quantities in 
the original observed map. Then we consider a pixelized flat 
2D image including CMB, Galactic foregrounds, instrumen- 
tal noise and extragalactic point sources; these four compo- 
nents are independent and are added up to form the image; 
whereas the CMB is Gaussian and the noise is considered 
as non-uniform but Gaussian-distributed at each pixel, the 
other two are non-Gaussian. 

Skewness and kurtosis can be obtained from the cumu- 
lants of the probability distribution function (pdf) p(x) . 
The characteristic function (j>{u) is defined as the Fourier 
transform of the pdf 



4>{u) = / p(x) e ux dx 
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and the cumulant function is K[u) = log(cf>(u)). The cumu- 
lants of order m, A m , can be calculated as derivatives of the 
cumulant function 



K m = 



1 d m K{u) 
i m du m 



(11) 



The second order cumulant is the variance, the third or- 
der cumulant is the third central moment {13 and the fourth 
order cumulant is fi4 — 3a 4 . The advantage of using the cu- 
mulants is that when the different components of the data 
are independent the cumulant of the sum is the sum of the 
cumulants. So, the cumulant of any order of the sum of 
CMB, noise, point sources and other foregrounds is the sum 
of the cumulants of the components. 

In a typical CMB experiment, the image is filtered with 
a Gaussian-shaped instrumental beam, this means that the 
map is convolved with a Gaussian function. This convolution 
changes the skewness and kurtosis of the original map. After 



the convolution, the temperature of a pixel is a linear com- 
bination of the temperatures of the neighboring pixels in the 
unfiltered map. For point sources, which are not correlated 
(this is true for the frequencies, 70 and 100 GHz, consid- 
ered i n this work, as can be seen in iGonzalez-Nuevo et al] 
(2005)), the cumulants of this linear combination can then 
be written as a linear combination of the cumulants of the 
original map. If we are filtering with a function ip(x), the cu- 
mulant of order m, K m , of the filtered image can be written 
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where A m is the original cumulant of order m in the pix- 
elized image and A p the pixel area. From this formula, we 
can obtain the cumulants of the Gaussian filtered image A„ 
from those of the original one. 
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± i n 



K„ 
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where I is the pixel size and 7 the beam dispersion, see eq. 
(4). We are going to filter now the image with the first three 
members of the MHWF, the whole process amounts to filter- 
ing the original image with a convolution of the Gaussian in- 
strumental beam and the corresponding wavelets. We could 
then apply Eq.(12) with ifi(x) this convolution. Finally, we 
can relate the cumulants K% of the wavelet filtered image 
to the same quantities K m in the original pixelized image. 
We write the formula for a general cumulant and a general 
member of index n of the MHWF 



K w — K 

I\ T „ — -fir? 
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in this formula /3 = R/j, R is the wavelet scale (see eq. 5) 
and the coefficients C mn are calculated as 



(tp n (x)) m x dx 



(15) 



This last expression can be easily computed for the first 
members of the MHWF and for m=3,4. 

When we consider the skewness of the total 
wavelet filtered image (CMB+point sources+Galactic fore- 
grounds-!- noise), it can be obtained as the total third order 
cumulant divided by o^, with a w the total rms deviation. 
The contribution of point sources to the total skewness of 
the wavelet filtered map is 



A' 



(16) 



where A™ is the third order cumulant (due to point sources) 
of the wavelet filtered map. The contribution of point sources 
to the total skewness of the Gaussian filtered map is 



h' 9 
^3 



(17) 



where A| is the third order cumulant (due to point sources) 
of the Gaussian filtered map and a g the total rms deviation 
of the Gaussian filtered maps. Finally, by using Eqs. (13) 
and (14) and the definition of the amplification, A, we can 
relate s w and s g . We obtain the following formula for the 



MHW1 
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for the MHW2, we obtain a similar formula 
22 

and finally for the MHW3 
124 °(1 + /3 2 )A% 



6561 



(19) 



(20) 



We want to point out that these contributions depend on the 
cubed amplification and will be higher for the wavelets with 
the highest amplifications, this suggests that MHW2 and 
MHW3 can be better at detecting non-Gaussianity due to 
point sources than MHW1. Besides, we will choose optimal 
scales when filtering with the MHWF so that the contribu- 
tion to the skewness is maximum. 

We can also obtain similar relations between the con- 
tribution of point sources k w to the kurtosis of the wavelet 
filtered maps and the contribution of point sources k g to 
the kurtosis of the beam filtered maps. So, we have for the 
MHW1 



(L-3) = ^(l + /3 2 )A 4 (fc 9 -3) 
for the MHW2 

(k. 



3) = |^(l + /3 2 )A 4 (fc 9 -3) 



8192 

and for the MHW3 



3) 



(21) 
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262144 

In this case, the contribution of point sources to the total 
kurtosis in the wavelet filtered map is proportional to A 4 . 



3 NON-GAUSSIANITY OF SIMULATED CMB 
MAPS 

In this section, we are going to consider the study of the 
non-Gaussianity in realistic simulated CMB maps. We have 
carried out simulations of CMB flat 2D maps of 12.8 x 12.8 
square degrees, generated with the cosmologica l parame- 
ters of the standard model iSpergel et al] 12003). We have 
then added the relevant Galactic foregrounds (free-free, syn- 
chrotron and dust emission) , we have used the Planck Refe- 
rence Sky, prepared by the members of the Planck Working 
Group 2. Galactic foregrounds have been selected from 12 
different regions of the sky (at Galactic latitude 6 > 30°). 
Each one of these 12 Galactic simulat ions has been used ac- 
cordin g to its representativeness as in lLopez-Caniego et al. I 
(2006). 

In order to give the most realistic numbers of detected 
extragalactic point sources, we hav e adop ted the model re- 
cently presented bv lDe Zotti et all (|2005l) . This new model 
takes into account the new data coming from high fre- 
quency radio surveys of extragalactic sources and discusses 
all source populations contributing to the number counts 
at these frequencies. It has thus proven capable of giving a 
better fit than before to all the currently published source 
number counts and relate d statistics at v ^ 10 GHz com- 
ing fr om different surveys ijRicci et alJl2004l : IWaldram et al] 
2003). This model also fits well th e number counts obta ined 
from the three year WMAP data l|Hinshaw et al.ll2006h . 

Since we are especially interested in applying our new 
method to the maps which will be provided by the Planck 



mission in the near future, and given that we can be very 
confident in the input source model counts, we have consid- 
ered the specific conditions of two Planck channels, which 
operate at 70 and 100 GHz. 1 

We are going to consider that the noise is anisotropic. 
Anisotropic noise has been generat ed using the simula - 
tion pipeline for the Planck mission dReinecke et al.l f2006). 
that has been made available for the Planck commu- 
nity by the Planck LevelS team. The s immission code 
Ivan Leeuwen et alJ 120021: IChallinor et alJ [2002) has been 
used to simulate the satellite dynamics during a mission 
flight of 14 months. The satellite will observe the microwave 
sky from the Lagrangian point L2, where it will be spin- 
stabilised by rotating at 1 r.p.m. The rotation axis will 
progress as the Earth moves around the Sun, therefore 
Planck will cover the full sky in half a year. The angle be- 
tween its optical and rotation axes will be slightly less than 
90°, thus Planck will scan the sky along circles. In order to 
cover the full sky, its rotation axis will make slow excursions 
from the Ecliptic plane with an amplitude of several degrees. 

For this simulation we have assumed a cycloidal scan- 
ning strategy with a 7° slow variation in the Ecliptic collat- 
itude of the z-axis. This scanning strategy implies that the 
sky will be covered inhomogeneously. As mentioned before, 
the sky is covered twice in a year. However, we have sim- 
ulated a flight duration of 14 months (the nominal mission 
duration). In the additional 2 months a region that covers 
approximately 1/3 of the total area of the sky is observed 
one more time. Figure 1 shows the observation pattern, in 
Ecliptic coordinates, for the 70 GHz channel. The observa- 
tion pattern is obtained from the satellite dynamics taking 
into account the position of the different detectors in the fo- 
cal plane, so the observation pattern for 100 GHz is slightly 
different from the one in Figure 1. In both cases we can see 
a broad "ring" where the number of observations is higher 
than in the rest of the sky. This "ring" covers ~ 1/3 of the 
sky. Since the final noise level is inversely proportional to 
the square root of the number of observations, the regions 
inside that ring have lower noise levels than the rest of the 
sky. In order to generate anisotropic noise patches for our 
work we have taken two representative zones of the sky, one 
inside the "low noise" ring and other outside it. The noise 
rms per pixel is obtained from the number of observations 
inside the pixel and the detector noise equivalent tempera- 
ture, that is specified in the Planck detector database. The 
final average noise rms levels per pixel (pixel size = 3') in 
thermodynamic AT/T units are 3.51 x 10~ 5 (1.35 x 10 -5 ) 
at 70 (100) GHz for the high noise zone and 2.19 x 10" 5 
(8.68 x 10" 6 ) at 70 (100) GHz for the low noise zone (inside 
the ring). To preserve the proportion of high to low noise 
regions in the sky we have simulated two high noise patches 
for every low noise patch. 

We consider 1000 simulations of CMB plus anisotropic 
noise at 70 and 100 GHz, these simulations are Gaussian and 
will serve as comparison with the non-Gaussian ones, one 



1 The characteristics of the channels, relevant for our purposes, 
are: a) pixel sizes, 3' and 3' at 70 and 100 GHz, respectively; 
b) Full Width Half Maximum (FWHM) of the circular gaussian 
beams, 14' and 9.5', respectively. In all the cases, we have used 
the estimated instrument performance goals available at the web 
site: http:/ /astro. estec.esa.nl/Planck 
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third of the simulations correspond to the low noise zone 
and two thirds to the high noise zone. Then, we generate 
1000 simulations including CMB, anisotropic noise, Galac- 
tic foregrounds and point sources for each of the following 

flux limits: Uy, 0.9 Jy, 0.8 Jy, , 0.1 Jy, i.e. we remove 

all the sources above the given flux limit (the third order 
and fourth order cumulants of the different simulated com- 
ponents can be seen in Table 2). We filter all the simulations 
with the MHW1, MHW2 and MHW3 at the optimal scales 
(obtained by maximizing the amplification) and calculate 
their skewness and kurtosis. 

Then, taking into account the skewness, we perform the 
following hypothesis test: we consider the null hypothesis, 
Ho, the simulation is Gaussian, against the alternative hy- 
pothesis, Hi, it is non-Gaussian. We set a significance level 
q = 0.05, this means that we reject the null hypothesis when 
a simulation with point sources has a skewness higher than a 
value ss, which corresponds to that of the 95% of the Gaus- 
sian simulations. The test is one-sided to the right, since the 
simulations with point sources have positive skewness. We 
define the power of the test as 1 — 5, with 5 the probabil- 
ity of accepting the null hypothesis when it is false, i.e. the 
power is the proportion of non-Gaussian simulations with a 
skewness higher than that of 95% of the Gaussian ones. The 
higher the power, the more efficient the method in detecting 
non- Gaussianity. 

We perform the same test with the kurtosis, setting a 
kurtosis limit ks; it is also one-sided to the right, since the 
kurtosis of point sources is higher than 3. 

The power of these tests for the different flux limits are 
shown in Table 1, they are also plotted in Figure 2. At 70 
GHz the power of the skewness test with a flux limit of 0.4 Jy 
is 8% (MHW1), 9% (MHW2) and 10% (MHW3). The power 
is higher for the kurtosis test : 11% (MHW1), 13% (MHW2) 
and 11% (MHW3) at 0.4 Jy. Therefore, for the MHW2 the 
probability that a patch with this flux limit has a kurtosis 
higher than that of 95% of the Gaussian simulations is 0.13. 
In practice, we could consider 126 patches (half sky) with 
the size of the simulated ones, if they are Gaussian the prob- 
ability of finding more than 10 with a kurtosis higher than 
fc 5 is 0.05, if we had point sources up to 0.4 Jy, this prob- 
ability would be 0.95, (this can be calculated by using the 
binomial distribution with probabilities 0.05 and 0.13 in the 
Gaussian and non-Gaussian case, respectively). Therefore, 
the detection of non-Gaussianity would be possible for this 
flux limit. 

For the 100 GHz channel, the results of the skewness 
test can also be seen in Table 1, for instance the power at 
0.3 Jy is 8% (MHW1), 13% (MHW2) and 12% (MHW3). 
The power is higher for the kurtosis test: 9% (MHW1), 18% 
(MHW2), 22% (MHW3) at 0.3 Jy, making detection possible 
down to this flux. The kurtosis test allows the detection of 
non-Gaussianity at 100 GHz for lower fluxes than at 70 GHz. 

We would like to comment that point sources could 
be detected in these channels, accepting a 5% of spuri- 
ous detect ions, down to 0.47 Jy (70 GHz) and 0.25 Jy 
(100 GHz ULopez-Caniego et al. Il2006ft . Therefore, even if 
we were able to perform a perfect subtraction (or mask- 
ing) of the point sources above these detection fluxes, we 
would still clearly detect a non-Gaussian signal due to un- 
subtracted point sources in the 70 GHz channel. For 100 
GHz, at the flux cut of 0.2 Jy (i.e below the theoretical detec- 



tion limit), the highest found power is 11% (corresponding 
to the kurtosis of the MHW3). Thus, it is likely that a non- 
Gaussian signal due to residual point sources is still present 
when analysing a presumably clean CMB map. These re- 
sults show the importance of studying and characterising 
the non-Gaussianity due to residual point sources (or other 
foregrounds) in order to avoid its misidentification as intrin- 
sic non-Gaussianity. 

In general the MHW2 and MHW3 perform better than 
the MHW1, producing higher powers; the kurtosis test is 
more powerful than the skewness one and the tests are more 
powerful at 100 GHz than at 70 GHz, (see Table 1 and Figure 
2). 

In order to assess the influence of the Galactic fore- 
grounds, we have also carried out 1000 simulations without 
point sources but including the Galactic foregrounds. The 
power in this case is always below 6% (see the no sources case 
in Table 1), confirming that these foregrounds contribute 
very little to the non-Gaussianity, when we filter with these 
wavelets at the optimal scales. 

We have also performed similar tests by filtering the 
maps with the Daubechies 4 wavelet, this wavelet has been 
alread y used in the study of the non-Gaussianity of CMB 
maps dBarreiro and Hobsonll200lh . The power of the skew- 
ness and kurtosis tests with this wavelet is never higher than 
50%, even for fluxes as high as 1 Jy and is always much lower 
than the power obtained with the MHWF, as can be seen 
in Figure 2. We have considered skewness and kurtosis of 
different detail coefficients (horizontal, vertical, diagonal) at 
different resolution levels and plotted in Figure 2 the power 
obtained for the detail coefficients which produce the highest 
power. It is clear that, whereas the members of the MHWF 
are adapted to the detection of point sources and the non- 
Gaussian signal they produce, this is not the case for the 
Daubechies family. 



4 CHARACTERISING POINT SOURCES 
4.1 Point source cumulants 

In the previous section, we have performed Gaussianity tests 
in order to detect the non- Gaussian signal generated by 
point sources in CMB maps. The estimation of the level of 
this signal is also of great interest. Therefore, it is worthwhile 
calculating the third order cumulant, K3, and the fourth or- 
der cumulant, K4, produced by point sources. The problem 
is that we could only know directly the values of K3 and K± 
of CMB maps including CMB, noise, point sources and other 
foregrounds, but not of the point source maps. A possible 
solution would be to filter the maps with suitable wavelets so 
that the main contribution to the cumulants would be that 
of point sources. As we have seen before, the MHWF at the 
optimal scales could perform well in this respect. Therefore, 
our idea is to filter CMB maps with the MHWF and to try 
to extract the parameters K3 and K4. However, we should 
first relate the cumulants of a wavelet filtered map to the cu- 
mulants of the original map, this can be done by using eqs. 
(14) and (15). These equations for the MHW1, MHW2 and 
MHW3 give us the relation between the original cumulants 
due to point sources and the wavelet filtered cumulants. 
Our goal is to estimate K3 and K4 from simulations 
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including CMB, noise, point sources and other foregrounds 
at 70 and 100 GHz. We filter these simulations with the 
MHW1, MHW2 and MHW3 at their optimal scales. We in- 
clude two additional scales for each wavelet to improve the 
fit, one arcminute less than the optimal scale and one ar- 
cminute more than the optimal scale. So for each simulated 
map we calculate nine quantities c l , the third order cumu- 
lants for the three scales and the three wavelets involved, 
where we have subtracted the average contribution of the 
anisotropic noise. We perform a chi-squared fit to the third 
order cumulant of the wavelet filtered maps g % Kz, with g l the 
coefficients, see eqs. (14) and (15), which relate the wavelet 
filtered third order cumulant with the original one K3. The 
chi-squared function can be written 



X 



(24) 



with dj the covariance matrix of the wavelet filtered third 
order cumulants obtained from 1000 simulations including 
CMB plus noise. The simulations used are flat patches of 
12.8 x 12.8 deg 2 similar to the ones explained in section 3. 
Minimizing the chi-squared expression with respect to K3, 
we estimate the third order cumulant due to point sources 
of the original map. We carry out 126 simulations (half sky) 
including CMB, anisotropic noise, point sources up to a cer- 
tain flux limit (from 0.1 Jy to 1 Jy) and other foregrounds. 
The third order and fourth cumulants of the different com- 
ponents of these simulations can be seen in Table 2. Note 
that in this Table, we write the cumulants after filtering the 
maps with the Gaussian instrumental beam, except for the 
noise. 

We obtain the point source cumulant K3 from each map 
through the chi-squared minimization. Finally, we calculate 
an average Kz representative of half sky. We perform the 
same process for the fourth order cumulant K4. The chi- 
squared function is in this case 



1 (<f - h 1 K 4 )(d J - h 3 K 4 ) 



(25) 



Now Dij is the covariance matrix of the fourth order cumu- 
lants, <T the fourth order cumulants of the wavelet filtered 
maps (subtracting the average anisotropic noise contribu- 
tion) and h? the coefficients which relate the wavelet filtered 
fourth order cumulants and the original ones. Minimizing 
with respect to K±, we estimate its value for each map. As 
before, we calculate an average value from 126 simulations 
(half sky) . Finally, we calculate the average value of K3 and 
Ki from ten half skies and their corresponding rms devia- 
tions. We have also calculated these cumulants directly from 
simulations (unfiltered) of the same size, which only include 
point sources up to a certain flux limit to compare our pre- 
vious estimations with the real values. The estimated values 
with the MHWF and the values obtained from the point 
source only simulations, K^ s and K% a , are presented in Ta- 
ble 3 for different flux limits from 0.1 Jy to 1 Jy. We have 
also calculated the relative error of Kz and K4 with respect 
to K 1 " and K 1 " for each half sky. We write in Table 3 the 
average and the dispersion of this error obtained from ten 
half skies. We can see in this table that the average relative 
error of K3 and K4, with respect to K^ s and iff s is below 
12% down to 0.6 Jy (70 GHz) and 0.4 Jy (100 GHz), so that 



values close to the real ones can be obtained down to low 
fluxes. 

Finally, in Table 4, we write K% s and K\ a calculated by 
carrying out point source simu lations in which we u se the 
Toffolatti source counts model (Toffola tti et alj|l9 98) as in- 
put for our source counts. It is clear that our estimations, Kz 
and K4, (see Tables 3 and 4) would allow us to distinguish 
between the De Zotti and the Toffolatti model. 



4.2 Fit of the counts to a power law 

Differential number counts can be modeled as a function of 
the flux S by a power law 

dn 
dS 



= AS~ 



(26) 



Though the fit is not very good for all the fluxes 
l|De Zotti et al.ll200.^ . it can be suitable for a certain range 
of fluxes, for instance for fluxes from 0.1 Jy to 1 Jy. Being 
able to obtain the values of the parameters A, the amplitude, 
and a, the slope, can be very useful to distinguish between 
different source counts models and to predict the number 
counts in the corresponding flux ranges. 

The third order and fourth order cumulants of the 
source temperature distribution in a pixelized map with 
pixel area A v can be easily calculated from the differential 
number counts assuming that sources are Poissonian dis- 
tributed. This is a good approximat ion for the frequencies 
(70 a nd 100 GHz) we are considering iGonzalez-Nuevo et alJ 
2005). In this case, the third order cumulant can be written 



K 3 = A. 



-2 3/ \ 

■ 9 0) 



, c c z dn 
dSS dS 



(27) 



where Su m is the flux limit and g is the conversion factor 
from fluxes to temperatures. 



{he) 2 (sinh(x/2)) 2 



(28) 



and x = hv fhsT. 

The formula for the fourth order cumulant is similar 

rSu m 



-3 4/ N 

. 9 W 







, c c 4 dn 

dSS dS' 



(29) 



If we assume a power law for the differential number 
counts, we could obtain Kz and K4 just by introducing the 
expression AS~ a in the previous formulas 



K 3 = 



K 



(4 - a) A 2 

Ag 4 s?- m a 



(30) 



(31) 



(5 - o)A% 

Then, A and a could be obtained from these equations, as- 
suming that we know K3, K4, and Su m - 

4 - 5c 



A = 



1 - c 

(4 - a)K 3 A 2 p 

3c(4-«) 
9 °lim 



with 



gKzSu r , 



(32) 
(33) 

(34) 
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After estimating K3 and K4 for the 70 and 100 GHz 
channels by using the fit with the MHWF (see the previous 
subsection), we can obtain the estimated values of the am- 
plitude, A, and the slope, a for each half sky by applying 
eqs. (32), (33) and (34). If we consider ten half skies, we can 
give an average value and the rms deviation. 

The results can be seen in Table 5, where they are com- 
pared with those obtained from the point source only sim- 
ulations; in this last case we have used the values K ps and 
K ps , calculated previously from the point source maps. The 
results at 70 GHz for a 1 Jy flux limit are: a = 2.19 ± 0.43, 
A = 24.3 ± 4.3, whereas we obtain from the point source 
simulations: a pa = 2.32 ± 0.06, A ps = 22.1 ± 1.5. For this 
flux limit, the average relative error of a {A) with respect to 
a ps ^ A p 3 j ig below 1Q% (15%) The estimation is not bad, 

but the rms deviation is high. We can obtain reasonable es- 
timations down to 0.7 Jy. It is clear that small errors in the 
determination of the cumulants produce high errors in the 
amplitude and slope. 

The values obtained from the fit for a flux limit of 1 Jy 
at 100 GHz are: a = 2.26±0.19, A = 21.0±4.2, whereas from 
the point source simulations we estimate: a ps = 2.35 ±0.16, 
A ps — 20.0 ± 3.7. The average relative error of a (A)with 
respect to a ps (A ps ) is below 5% (6%). We obtain reasonable 
values down to 0.7 Jy, with less dispersion than for the 70 
GHz channel. 

In Figure 3, the fit for both channels at 1 Jy is compared 
to the differential number counts of the De Zotti model, 
which was used for the simulations and with the estimation 
from the point source simulations. As can be seen in the 
plots, the fit is quite good for fluxes between 0.1 Jy and 1 

Jy. 



5 CONCLUSIONS 

We have considered a natural generalization of the MHW1 
on the plane, R 2 , based on the iterative application of the 
Laplacian operator to the MHW1. We have called this group 
of wavelets the Mexican H a t Wavelet Family (MHWF ) 
l|Gonzalez-Nuevo et alJ hoOrl lLopez-Caniego et al. I EoOrl) . 
We have applied the MHWF to the study of the non- 
Gaussianity produced in CMB maps by extragalactic point 
sources, a very important issue for the separation and anal- 
ysis of different non-Gaussian components in the CMB. 

We have first derived a formula, eq. (14), which re- 
lates the cumulants of 2D point source maps filtered with 
the MHWF with the cumulants of the original point source 
maps. From this formula, we have obtained several equa- 
tions, eqs. (18-23), which relate the contribution of point 
sources to the skewness and kurtosis of wavelet (MHWF) 
filtered 2D CMB maps to their contribution to the skewness 
and kurtosis of beam filtered CMB maps. 

According to these formulas, the skewness of a wavelet 
filtered map is proportional to A 3 , with A the amplification, 
and the kurtosis is proportional to A 4 . This means that the 
higher the amplification the higher the non-Gaussian signal 
will be. Therefore, the MHW2 and MHW3 at their optimal 
scales can be very suitable filters to detect the non-Gaussian 
signal due to point sources, s ince they produce a higher 
amplification than the MHW 1 l)Gonzalez-Nuevo et al Jl2006t : 
lLopez-Caniego et al. 1 l2006l) . 



We have then considered Gaussian simulations of CMB 
plus anisotropic noise and compared them with simulations 
including all the relevant components: CMB, anisotropic 
noise, point sources and Galactic foregrounds. The simula- 
tions are 2D maps of 12.8 x 12.8 square degrees with the char- 
acteristics of the Planck experiment at 70 GHz and 100 GHz. 
We perform a Gaussianity test based on the skewness and 
kurtosis of the maps filtered with the MHW1, MHW2 and 
MHW3 at the optimal scales (those which produce the high- 
est amplifications), we compare 1000 simulations of CMB 
plus noise with 1000 simulations including point sources up 
to different flux limits from 0.1 Jy to 1 Jy. 

The power of the test is higher for the MHW2 and the 
MHW3 than for the standard MHW1, the kurtosis test is 
more powerful than the skewness test. We could detect the 
non- Gaussianity due to point sources down to a flux limit of 
0.4 Jy (70 GHz) and 0.3 Jy (100 GHz), the powers are higher 
for 100 GHz than for 70 GHz. The results can be seen in Fi- 
gure 2 and Table 1. We woul d like t o point out that accord- 
ing to lLopez-Caniego et al. I i2006T) point sources could be 
detected in these channels, assuming a 5% of spurious detec- 
tions, down to 0.47 Jy (70 GHz)and 0.25 Jy (100 GHz). This 
means that, using the MHWF, we can see a non-Gaussian 
signal due to weak point sources, which are in the limit of 
detection of current source extraction techniques. Therefore, 
it is very important to characterise the non-Gaussianity due 
to unsubtracted point sources, in order to avoid its possible 
misidentification as intrinsic non-Gaussianity. 

The influence of Galactic foregrounds is negligible, since 
when we consider simulations without point sources, but 
keeping the Galactic foregrounds, the powers go down below 
6%. We have also studied the performance of the Daubechies 
4 wavelet, obtaining powers below 50% even for flux limits 
as high as 1 Jy. The MHWF performs much better than this 
wavelet (see Figure 2), since it is adapted to deal with point 
sources. 

This method could be used for all the frequencies of the 
Planck experiment and could also be extended to spherical 
CMB images, as the ones provided by the Planck experi- 
ment. 

We have also tried to estimate the third order cumulant 
Kz and the fourth order cumulant K4, produced by point 
sources at the 70 and 100 GHz Planck channels. We filter 
12.8 x 12.8 square degrees CMB maps with the characteris- 
tics of the Planck experiment with the MHW1, MHW2 and 
MHW3 at three scales for each wavelet (including the opti- 
mal ones). The maps include CMB, anisotropic noise, point 
sources and other foregrounds. We calculate the third and 
fourth order cumulants of these wavelet filtered maps and 
write chi-squared functions, eqs. (24) and (25), which allow 
us to estimate the cumulants of the original point source 
maps; we have used eq. (14) to relate the cumulants of the 
wavelet filtered maps with those of the original maps. We 
have compared our estimated values with those, K ps and 
K ps , obtained directly from simulations of the same type 
but including only unfiltered point sources without any other 
component. The results can be seen in Table 3. The aver- 
age error (from ten half skies) is below 12% down to a flux 
limit of 0.6 Jy (70 GHz) and 0.4 Jy (100 GHz). The deter- 
mination of these cumulants would allow us to distinguish 
between different source counts models (see Tables 3 and 4) . 

Finally, we have also studied the fit of the differen- 
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tial number counts at 70 and 100 GHz to a power law 
dn/dS = AS~ a . From the cumulants previously estimated, 
we calculate by using eqs. (32-34) the amplitude A and slope 
a of the differential number counts. The results at 70 GHz 
for a 1 Jy flux limit are: a = 2.19 ± 0.43, A = 24.3 ± 4.3. 
When we calculate these magnitudes from the cumulants 
K^ s and K ps , computed directly from the point source only 
simulations, we obtain a ps = 2.32 ± 0.06, A ps = 22.1 ± 1.5. 
These are the average and rms deviation from ten half skies. 
For this flux limit, the average relative error of a (A) with 
respect to a ps {A pa ) is below 10% (15%). The dispersion in 
the values obtained with the MHWF fit is high, though the 
average is consistent with the value obtained directly from 
point sources down to a flux limit of 0.7 Jy, as can be seen 
in Table 5. 

When we consider the 100 GHz channel, the values 
obtained from the MHWF fit for a flux limit of 1 Jy are: 
a = 2.26 ± 0.19, A = 21.0 ± 4.2, the results from the point 
source simulations are: a ps = 2.35 ± 0.16, A pa = 20.0 ± 3.7. 
The fit is better for this channel than for the 70 GHz one, 
having a lower dispersion. For this flux limit, the average 
relative error of a (A) with respect to a pa (A ps ) is below 
5% (6%). The results are consistent down to a flux limit of 
0.7 Jy. The detailed values of the parameters a and A can 
be seen in Table 5. 

In Figure 3 we plot the differential source counts of the 
De Zotti model, the fit to a power law by using the MHWF 
and the fit obtained from the point source only simulations, 
in both cases with a 1 Jy flux limit. As we can see in this 
figure, both fits are consistent with the input counts in the 
range between 0.1 Jy and 1 Jy. 
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Table 1. Power (percentage) of the skewness and kurtosis tests at 70 and 100 GHz for different flux limits, for the first three members 
of the MHWF. In the first, second and third columns, we show the power of the skewness test for the MHW1, MHW2 and MHW3, 
respectively. In the fourth, fifth and sixth columns, we write the power of the kurtosis test for the MHW1, MHW2 and MHW3. 
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Table 2. Third order cumulant K3 and fourth order cumulant K4 of the different components of our simulations: noise (high noise 
zone), noise (low noise zone), CMB, Galactic foregrounds and point sources, up to flux limits of 1 Jy, 0.5 Jy and 0.1 Jy. The cumulants 
are obtained from the simulations filtered with the Gaussian beam, except for the noise. We show the average and rms deviation from 
ten half skies. 
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(2.38+0.26) xl0~ 


16 


(2.37+0.07) xl0~ 


16 


7.8+5.8 


0.5 Jy 


(1. 13+0. 24) x 10" 


13 


(1.12±0.04)xl0~ 


13 


18.0+10.0 


(1.14+0.42) x 10" 


16 


(1.24+0.06) x 10" 


16 


26.0+23.1 


0.4 Jy 


(8. 37+2. 85) x 10" 


14 


(8.80+0.42) xlO~ 


■14 


25.2+18.1 


(9.54+3.03) xl0~ 


17 


(8.56+0.56) x 10" 


■17 


26.7+23.2 


0.3 Jy 


(6. 44+2. 24) x 10" 


14 


(4.33+0.19) xlO~ 


■14 


58.4+46.7 


(4.26+2.58) x 10" 


17 


(3.03+0.16) x 10" 


■17 


77.6+61.7 


0.2 Jy 


(1.74+3.15) x 10" 


14 


(2.08+0.08) xlO~ 


■14 


120+86 


(9.89+17.71) x 10" 


-18 


(1.02+0.05) x 10" 


■17 


134+110 


0.1 Jy 


(1.09+3. 11) x 10" 


14 


(3.90+0.10) xlO~ 


■15 


586+538 


(7.90+34.6) x 10" 


18 


(8.56+0.29) xlO" 


19 


3350+2037 


100 GHz 


K 3 




K T 




Aft 


K A 




Kl s 




AK 4 


1 Tv 


(6 90+0 54")xl0 _ 


■14 


(6 98+0 40"lxl0~ 


■14 


2.3+1.7 


C8 45+0 76 1 ) X 1 0" 


17 


(8 38+0 60~)x10" 


■17 


2.5+2.3 


0.9 Jy 


(6.39+0.50) x 10~ 


14 


(6.54+0.41) x 10" 


■14 


2.8+2.4 


(7.42+0.69) x 10" 


17 


(7.50+0.64) x 10" 


■17 


2.5+1.6 


0.8 Jy 


(5.38+0.59) x 10~ 


14 


(5.50+0.42) x 10" 


■14 


4.6+3.2 


(5.90+0.65) x 10" 


17 


(5.83+0.58) x 10" 


■17 


3.6+2.3 


0.7 Jy 


(4.29+0.32) x 10~ 


14 


(4.25+0.24) x 10" 


■14 


4.3+2.7 


(3.86+0.29) x 10" 


17 


(3.88+0.29) x 10" 


■17 


1.9+1.5 


0.6 Jy 


(2.80+0.14) x 10" 


14 


(2.99+0.17) x 10" 


■14 


6.0+5.3 


(2.36+0.16) x 10" 


17 


(2.31+0.17) x 10" 


■17 


4.6+3.4 


0.5 Jy 


(1.75+0.13) x 10" 


14 


(1.88+0.09) x 10" 


■14 


7.2+7.8 


(1.23+0.07) x 10" 


17 


(1.17+0.07) x 10" 


■17 


6.1+3.8 


0.4 Jy 


(1.31+0.20) x 10" 


14 


(1.43+0.10) x 10" 


■14 


11.3+7.4 


(7.97+0.41) x 10" 


18 


(7.79+0. 59) x 10" 


■18 


4.0+2.7 


0.3 Jy 


(6.17+1.17) x 10" 


15 


(7.19+0.27) x 10" 


■15 


16.6+14.9 


(3.35+0.39) x 10" 


18 


(2.81+0.13) x 10" 


■18 


22.8+11.1 


0.2 Jy 


(1.97+1.96) x 10" 


■15 


(3.39+0.09) x 10" 


■15 


57.2+40.4 


(1.52+0.25) x 10" 


18 


(0.93+0.03) x 10" 


■18 


63.4+27.2 


0.1 Jy 


(-8.30+15.4) x 10" 




(6.10+0.16) x 10" 


■16 


282+187 


(5.43+3.85) x 10" 


19 


(7.57+0.26) x 10" 


■20 


622+476 



Table 3. First and fourth columns: third order cumulant K3 and fourth order cumulant K4 estimated from the CMB maps by carrying 
out a fit with the MHWF, the average from ten half skies and the rms deviation are written. Second and fifth columns: third order 
cumulant K ps and fourth order cumulant K PS estimated from the the point source simulations (average and rms deviation from ten half 
skies). Third and sixth columns: relative error (percentage) of K3 and K4 with respect to K ps and K ps (average and rms deviation from 
ten half skies). 



70 GHz 






Kf 




1 Jy 


(8.6+0.4) x 10" 


13 


(1.8+0.1) x 10" 


15 


0.9 Jy 


(8.1+0.3) x 10" 


13 


(1.6±0.1)xl0" 


15 


0.8 Jy 


(7.1±0.2)xl0" 


13 


(1.32+0.06) x 10 


-15 


0.7 Jy 


(5.4+0.2) x 10" 


13 


(8.7±0.5)xl0" 


Ifi 


0.6 Jy 


(3.9±0.1)xl0" 


13 


(5.2±0.2)xl0" 


Ifi 


0.5 Jy 


(2.60+0.06) x 10 


-13 


(2.8±0.1)xl0" 


Ifi 


0.4 Jy 


(2.02+0.05) x 10 


-13 


(1.90+0.07) x 10 


-16 


0.3 Jy 


(1.07+0.02) x 10 


-13 


(7.3±0.2)xl0" 


17 


0.2 Jy 


(5.4±0.1)xl0" 


14 


(2.59+0.06) x 10 


-17 


0.1 Jy 


(1.17+0.01) x 10 


-14 


(2.55+0.03) x 10 


-18 


100 GHz 


K ps 
"■3 




K p3 




1 Jy 


(1.49+0.07) x 10 


-13 


(1. 7+0.1) x 10~ 


16 


0.9 Jy 


(1.4+0.1) x 10" 


13 


(1.6±0.1)xl0- 


Ifi 


0.8 Jy 


(1.20+0.04) x 10 


-13 


(1.24+0.06) x 10 


-16 


0.7 Jy 


(9.1±0.3)xl0" 


14 


(8.1±0.4)xl0" 


17 


0.6 Jy 


(6.6±0.2)xl0" 


14 


(5.0±0.2)xl0" 


17 


0.5 Jy 


(4.3±0.1)xl0- 


14 


(2.6±0.1)xl0" 


17 


0.4 Jy 


(3.4±0.1)xl0" 


14 


(1.81+0.05) x 10 


-17 


0.3 Jy 


(1.81+0.03) x 10 


-14 


(6.8±0.1)xl0" 


18 


0.2 Jy 


(9.3±0.1)xl0" 


15 


(2.46+0.04) x 10 


-18 


0.1 Jy 


(1.98+0.02) x 10 


-15 


(2.41+0.03) x 10 


-19 



Table 4. Third order cumulant K ps and fourth order cumulant K ps calculated from the the point source simulations carried out with 
the source number counts of the Toffolatti et al. model (average and rms deviation from ten half skies). 
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70 GHz 


a 




Act 


A 


A ps 


AA 


i Jy 

0.9 Jy 
0.8 Jy 
0.7 Jy 


2.19±0.43 
2.08±0.50 
1.81±0.72 
2.24±0.46 


2.32±0.06 
2.10±0.12 
2.03±0.12 
2.06±0.05 


9.9±16.1 
15.5±13.4 
27.2±21.1 
18.6±12.3 


24.3±4.3 
28.3±9.4 
28.3±9.4 
25.2±9.4 


22.1±1.5 
27.0±3.6 
27.0±3.6 
28.8±2.6 


14.9±14.9 
18.9±17.2 
44.6±43.3 
26.5±15 


100 GHz 


a 




Act 


A 




AA 


i Jy 

0.9 Jy 
0.8 Jy 
0.7 Jy 


2.26±0.19 
2.16±0.14 
1.84±0.26 
1.93±0.16 


2.35±0.16 
2.23±0.13 
2.07±0.12 
2.14±0.11 


4.2±3.5 
4.3±3.8 
10.9±13.7 
10.0±6.6 


21.0±4.2 
23.4±3.4 
30.6±5.7 
28.6±3.7 


20.0±3.7 
22.9±3.2 
26.5±4.1 
24.8±3.8 


5.6±5.9 
5.7±4.9 
17.5±22.1 
16.4±12.6 



Table 5. First and fourth columns: slope a and amplitude A of the differential number counts estimated from the CMB maps by 
carrying out a fit with the MHWF. Second and fifth columns: slope a pa and amplitude A ps estimated from the point source simulations. 
Third and sixth columns: relative error (percentage) of a and A with respect to a ps and A ps . In all the cases, we write the average and 
rms deviation from ten half skies 
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Figure 1. Observation pattern for the 70 GHz channel. The map 
shows the logarithm of the number of hits (number of times the 
detector observes a given pixel during the flight) in Ecliptic co- 
ordinates. 



14 Argiieso et al. 





Figure 2. Power against the flux limit for the skewness test at 70 GHz (top left panel), the kurtosis test at 70 GHz (top right panel), 
the skewness test at 100 GHz (bottom left panel) and the kurtosis test at 100 GHz (bottom right panel). The solid line represents the 
MHW1 results, the dashed line the MHW2 results, the dotted line the MHW3 results and the dash-dotted line the Daubechies4 results. 
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S(Jy) 




S(Jy) 



Figure 3. Top panel: (70 GHz channel) differential number counts against the flux in Jy according to the De Zotti model (solid line), 
the fit to a power law by using the MHWF (dashed line) with the average value of the parameters for a flux limit of 1 Jy and the same 
fit by using the original point source maps (dotted line) with the average value of the parameters for a flux limit of 1 Jy. Bottom panel: 
the same for the 100 GHz channel. 



